Ejercicio 1

Para nuestro trabajo hemos elegido la base de datos pública, que hemos encontrado en el sitio web kaggle.com:

https://www.kaggle.com/uciml/biomechanical-features-of-orthopedic-patients (data recodida inicialmente de Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science)

Hemos eligido este data set por: - Una de las problemas frecuentes en la medicina en el ambito contemporaneo es la optimización de recursos donde sea posible. Y una área de aplicación de la analísis es en la etápa del diagnóstico de algunas enfermedades en las etapas tempranas, así por ejemplo, el diagnóstico de la hernia intervertebral exige utilizar la prubas caras y que tardan mucho en realizarse, cuando hay alternativas, que, quizas, no son tan fiables, pero puedan ayudar a la hora de selecionar pacientes para la siguiente etápa de estudios. En este caso se puede hablar de comparación entre estudio con tomografia computarizada vs la resonancia magnética. - Es un data-set público y facil de acceder. - Es un data-set relativamente límpio y ya preparado para realizar analisis de datos sobre el.

Ejercicio 2

Esta dataset tiene 2 archivos .csv (column_2C_weka y column_3C_weka), que se puede juntar y valorar conjuntamente: Es una base de datos que contiene informacaión sobre 310 personas, en las que evaluaron posición de vertebras de la columna lumbar utilizando varios parametros (tabla 1) y exolicación visual con la imagen del articulo^1. Dentro de esta muestra habia 100 personas que consideramos “sanos” o “normales” y 210 (“anormales”) que consideramos que tenian algun tipo de patología: “hernia” o “spondilolistesis”.

tabla 1:

Variable Expicación del variable
pelvic_incidence, numerico en grados Angulo entre la línea que pasa del centro de la lámina terminal del S1 al centro de las cabezas femorales, y la línea perpendicula a la lámina terminal de S1
pelvic_tilt, numerico en grados Angulo entre la línea que pasa del centro de la lámina terminal del S1 al centro de las cabezas femorales y una línea vertical
lumbar_lordosis_angle, numerico angulo entre la lámina terminal superior del L1 y S1
sacral_sclope, numerico en grados angulo entre la línea perpendicula a la lámina terminal de S1 y una línea horizontal
pelvic_radius, numerico en grados Es un angulo entre la línea vertical y una línea del centro de las cabezas femorales hasta el borde posterior de la lámina terminal de S1
degree_spondylolysthesis, numerico en mm distancia de desplazamiento de la vertebra superior sobre la vertebra inferior (positivo - anterior y negativo - posterior)
class es posible utilizar este variable como character o factor.
class_2c, factor Classificación de data como “normal” y “anoraml”
class_3c, factor Classificación de data como “normal”, “Hernia” y “spondilolistesis”

En esta imagen presentaremos la visualización de lo angulos medidas en la base de datos. Lumbar_angles

Las bibliografía:

  1. Sergides IG, McCombe PF, White G, Mokhtar S, Sears WR. Lumbo-pelvic lordosis and the pelvic radius technique in the assessment of spinal sagittal balance: strengths and caveats. Eur Spine J. 2011;20 Suppl 5(Suppl 5):591-601. doi:10.1007/s00586-011-1926-z

  2. https://www.orthobullets.com/spine/2071/lumbar-spine-anatomy para la fecha de 25.12.2021

# Vamos a introducir los datos en RStudio. Para esto vamos a utilizar la libreria básica del R 
column_2C <- read_csv("archive/column_2C_weka.csv", 
                    col_types = cols(pelvic_incidence = col_number(), 
                                           `pelvic_tilt numeric` = col_number(), 
                                           lumbar_lordosis_angle = col_number(), 
                                           sacral_slope = col_number(), pelvic_radius = col_number(), 
                                           degree_spondylolisthesis = col_number()))
column_3C <- read_csv("archive/column_3C_weka.csv", 
                      col_types = cols(pelvic_incidence = col_number(), 
                                             pelvic_tilt = col_number(), lumbar_lordosis_angle = col_number(), 
                                             sacral_slope = col_number(), pelvic_radius = col_number(), 
                                             degree_spondylolisthesis = col_number()))
# Vamos a ajustar los nombres de variables y confirmamos que son iguales
colnames(column_2C) <- colnames(column_3C)
colnames(column_2C) %>% 
  set_names() %>% 
  map(~(column_2C[,.x]==column_3C[,.x])) %>% 
  as.data.frame() %>%
  summary() #utilizaremos función map de la libreria purr de tidyverse, aplicando la función a todos los elementos de las columnas correspondientes.
##  pelvic_incidence pelvic_tilt    lumbar_lordosis_angle sacral_slope  
##  Mode:logical     Mode:logical   Mode:logical          Mode:logical  
##  TRUE:310         TRUE:310       TRUE:310              TRUE:310      
##                                                                      
##  pelvic_radius  degree_spondylolisthesis   class        
##  Mode:logical   Mode:logical             Mode :logical  
##  TRUE:310       TRUE:310                 FALSE:210      
##                                          TRUE :100
# Y cambiamos el nombre de la columna "class" para poder juntar los dos data frame.
column_2C<-rename(column_2C, "class_2" = "class")
column_3C<-rename(column_3C, "class_3" = "class")

# Creamos un data frame con información de los dos archivos iniciales.
ortho_df <- merge(column_2C, column_3C, by = c("pelvic_incidence", "pelvic_tilt", "lumbar_lordosis_angle", "sacral_slope", "pelvic_radius", "degree_spondylolisthesis"))
# Convertimos la varibel class en un factor, con los niveles que van de nivel normal a nivel con alteraciones más faciles de encontrar (no se puede evidenciar la hernía intervertebral en la radiografia, pero la espondilolistesis - sí).
ortho_df <- ortho_df %>%
  mutate("class_2" = factor(ortho_df$`class_2`,
                            levels = c("Normal", "Abnormal")))
ortho_df <- ortho_df %>%
  mutate("class_3" = factor(ortho_df$`class_3`,
                            levels = c("Normal", "Hernia", "Spondylolisthesis")))
# En algunas ocasiones, para el objetivo de visualizar la data con mayor facilidad, vamos a cerar un data set con 
ortho_df2 <- ortho_df %>%
  pivot_longer(cols = c(pelvic_incidence, pelvic_tilt, pelvic_radius, sacral_slope, lumbar_lordosis_angle),
               names_to = "type_of_measure", values_to = "angles") %>% 
  select(type_of_measure, angles, everything())
ortho_df3 <- ortho_df %>% filter(degree_spondylolisthesis<400) # limpiamos el data.frame de un valor evidentemente erroneo (degrre_spondilolistesis = 418 mm)

Ejercicio 3:

De todas las cuestiones que podemos realizar con estos datos, seleccionaremos las dos que, para nosotros, son las más destacables y extrapolables a una utilidad médica real.

La primera cuestión sería relativa a la comparación de grupos. Pregunta 1: ¿Existen diferencias significativas entre las variables angulares del grupo “Normal” respecto al grupo “Abnormal” (personas con hernia o con espondilolistesis)?. En otras palabras, ¿el hecho de tener hernia o espondilolistesis afecta significativamente en alguna de las variables de estudio? Su respuesta nos permitirá discernir qué métodos de medida difieren estadísticamente respecto al grupo “Normal” y llegar a un mejor acercamiento a nivel de diseño de prótesis o clasificación de características de pacientes en bases de datos.

Para resolver esta pregunta compararemos los datos de los grupos a lo largo de este estudio, visualizaremos los estadísticos básicos de nuestra muestra y realizaremos un análisis de varianza (ANOVA) de los grupos (véase los ejercicios 4 y 7).

Pregunta 2: ¿Existe alguna relación o correlación entre las variables del estudio? Esta pregunta nos permitirá conocer la relación entre las variables de estudio para saber si puede existir cierta redundancia en la medida de la anatomía de la columna lumbar. Además, nos permitirá saber si existe un patrón significativo de influencia de una variable determinada a otra/otras variable/s.

Esta cuestión podrá resolverse mediante los estudios de regresión líneal y regresión múltiple (véase los ejercicios 5 y 6). Además, se realizará un análisis de clústering (véase ejercicio 7) para observar si es posible agrupar los datos de nuestra muestra.

Ejercicio 4:

El primer elemento de nuestros datos que analizaremos es la estadística básica y el comportamiento de cada variable entre los grupos.

Primeramente observaremos la estadística básica de cada variable mediante el comando “summary()”. Éste nos muestra los valores mínimos, máximos, primer y tercer cuartil, mediana y media para cada variable. Más adelante calcularemos su desviación estándar para cada grupo.

Al ser una variable de tipo caracter, mostramos las observaciones por grupos en las variables “class_2” y “class_3” mediante la herramienta “table()”.

summary(ortho_df)
##  pelvic_incidence  pelvic_tilt     lumbar_lordosis_angle  sacral_slope   
##  Min.   : 26.15   Min.   :-6.555   Min.   : 14.00        Min.   : 13.37  
##  1st Qu.: 46.43   1st Qu.:10.667   1st Qu.: 37.00        1st Qu.: 33.35  
##  Median : 58.69   Median :16.358   Median : 49.56        Median : 42.40  
##  Mean   : 60.50   Mean   :17.543   Mean   : 51.93        Mean   : 42.95  
##  3rd Qu.: 72.88   3rd Qu.:22.120   3rd Qu.: 63.00        3rd Qu.: 52.70  
##  Max.   :129.83   Max.   :49.432   Max.   :125.74        Max.   :121.43  
##  pelvic_radius    degree_spondylolisthesis     class_2   
##  Min.   : 70.08   Min.   :-11.058          Normal  :100  
##  1st Qu.:110.71   1st Qu.:  1.604          Abnormal:210  
##  Median :118.27   Median : 11.768                        
##  Mean   :117.92   Mean   : 26.297                        
##  3rd Qu.:125.47   3rd Qu.: 41.287                        
##  Max.   :163.07   Max.   :418.543                        
##               class_3   
##  Normal           :100  
##  Hernia           : 60  
##  Spondylolisthesis:150  
##                         
##                         
## 
table(ortho_df$class_2)
## 
##   Normal Abnormal 
##      100      210
table(ortho_df$class_3)
## 
##            Normal            Hernia Spondylolisthesis 
##               100                60               150

Obsérvese que podemos determinar los límites anotados para cada variable y el comportamiento general mediante la media. Estos valores nos servirán más adelante para determinar un intervalo para presentar las variables en gráficos.

La variable “pelvic_incidence” tiene valores de 26.15º hasta 129.83º.

La variable “pelvic_tilt” tiene valores de -6.555º hasta 49.432º.

La variable “lumbar_lordosis_angle” tiene valores de 14º hasta 125.74º.

La variable “sacral_slope” tiene valores de 13.37º hasta 121.43º.

La variable “pelvic_radius” tiene valores desde 70.08º hasta 163.07º.

La variable “degree_spondylolisthesis” tiene valores desde -11.058 mm hasta 418.543 mm.

Las variables “class_2 y class_3” incluyen 210 personas en el grupo “Abnormal” (60 con hernia y 150 con espondilolistesis) y 100 personas en el grupo “Normal”.

Aunque esta descripción nos sirve para acotar los valores de las variables con los que trabajamos, los datos comprenden todos grupos de estudio en conjunto. Nuestro objetivo principal es la comparación entre personas del grupo “Normal” (sin anomalías aparentes relativas a la postura pélvica y lumbar) y “Abnormal” (personas con hernia o espondilolistesis diagnosticada) para determinar si existen diferencias significativas en las variables.

El hecho de dividir los datos podría llegar a mostrar una diferencia en el comportamiento de los mismos. Es por eso que realizaremos un análisis de comparación de variables según los grupos. El valor que compararemos será la media de cada variable, con su desviación estándar anotada. Asimismo, para hacer un análisis visual, las diferencias entre grupos de estudio se representarán en gráficos de cajas (boxplots). En cada caso se utilizarán los intervalos calculados anteriormente para que la representación de cada gráfico esté en el mismo intervalo y que así su comparación sea más intuitiva y fácil.

Hace falta remarcar que este análisis lo haríamos para las variables de tipo angular (las cinco primeras variables). No obstante, al tratarse de un trabajo de ámbito educativo, para la variable “grado de espondilolistesis” realizaremos la comparación igualmente, pero de antemano sabemos que será superior en pacientes del grupo “espondilolistesis” respecto al resto, ya que tienen dicha afectación diagnosticada.

Mediante la herramienta “aggregate” calculamos la media y la desviación estándar para cada grupo de nuestras variables.

aggregate(ortho_df[, 1:6], list(ortho_df$class_3), mean)
##             Group.1 pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## 1            Normal         51.68524    12.82141              43.54260
## 2            Hernia         47.63841    17.39879              35.46352
## 3 Spondylolisthesis         71.51422    20.74804              64.11011
##   sacral_slope pelvic_radius degree_spondylolisthesis
## 1     38.86383      123.8908                 2.186572
## 2     30.23961      116.4750                 2.480251
## 3     50.76619      114.5188                51.896687
aggregate(ortho_df[, 1:6], list(ortho_df$class_3), sd)
##             Group.1 pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## 1            Normal         12.36816    6.778503             12.361388
## 2            Hernia         10.69713    7.016708              9.767795
## 3 Spondylolisthesis         15.10934   11.506169             16.397068
##   sacral_slope pelvic_radius degree_spondylolisthesis
## 1     9.624004      9.014246                 6.307483
## 2     7.555388      9.355720                 5.531177
## 3    12.318813     15.579995                40.108030

Representaremos estos estadísticos a continuación para cada variable.

Variable 1: Incidencia pélvica (pelvic_incidence)

par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Incidencia pélvica", xlab="Normal", ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_2=="Abnormal"],col="Tomato",main="Incidencia pélvica", xlab="Abnormal",ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Hernia"],col="coral1",main="Incidencia pélvica", xlab="Hernia", ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Incidencia pélvica", xlab="Spondylolisthesis", ylim=c(20,130))

par(mfrow=c(1,1))

Observamos que la media del grupo “Normal” es ligeramente inferior a la del grupo “Abnormal”. Al haber más pacientes en la muestra con “Spondylolisthesis” que en la de “Hernia”, la media del grupo “Abnormal” resulta ser superior a la del grupo “Normal”, con una desviación estándar muy alta. Este hecho se observa en los gráficos, donde parece que el grupo “Spondylolisthesis” supera significativamente la media otorgada por el grupo “Normal”. Como sólo observamos tres outliers aparentes en el grupo de espondilolistesis (la muestra no parece tener mucha distorsión), nuestra primera hipótesis será que los pacientes con espondilolistesis podrían tener una mayor incidencia pélvica general. Así, es probable que la variable “incidencia pélvica” tenga cierto peso a la hora de analizarse en dicho grupo.

Variable 2: Inclinación pélvica (pelvic_tilt)

par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Inclinación pélvica", xlab="Normal", ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_2=="Abnormal"],col="Tomato",main="Inclinación pélvica", xlab="Abnormal",ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Hernia"],col="coral1",main="Inclinación pélvica", xlab="Hernia", ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Inclinación pélvica", xlab="Spondylolisthesis", ylim=c(-7,50))

par(mfrow=c(1,1))

Se observa que la media del grupo “Normal” para la inclinación pélvica es ligeramente inferior que la de los otros grupos. No obstante, al existir una desviación estándar realmente alta para todos los grupos en esta variable (es decir, existe una gran variabilidad), es probable que la distinción entre grupos no llegue a ser significativa. Por lo tanto, no esperamos diferencias significativas entre individuos de los distintos grupos para esta variable.

Variable 3: Ángulo de lordosis lumbar (lumbar_lordosis_angle)

par(mfrow=c(1,4))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Lordosis lumbar", xlab="Normal", ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_2=="Abnormal"],col="Tomato",main="Lordosis lumbar", xlab="Abnormal",ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Hernia"],col="coral1",main="Lordosis lumbar", xlab="Hernia", ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Lordosis lumbar", xlab="Spondylolisthesis", ylim=c(13,126))

par(mfrow=c(1,1))

En este caso observamos que la media del grupo “Normal” es parecida a la del grupo “Hernia”, pero difiere bastante de la del grupo “Spondylolisthesis”. Aunque exista algún outlier apartado que pueda hacer aumentar la desviación estándar de cada grupo (véase los puntos superiores que se distancian de la media en cada gráfico), es probable que esta comparativa pueda ser significativa. Así, nuestra segunda hipótesis es que los pacientes con espondilolistesis podrían tener un mayor ángulo de lordosis lumbar. Y por lo tanto, la variable “ángulo de lordosis lumbar” podría tener más peso en este grupo de pacientes.

Variable 4: Pendiente sacral (sacral_slope)

par(mfrow=c(1,4))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Pendiente sacral", xlab="Normal", ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_2=="Abnormal"],col="Tomato",main="Pendiente sacral", xlab="Abnormal",ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Hernia"],col="coral1",main="Pendiente sacral", xlab="Hernia", ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Pendiente sacral", xlab="Spondylolisthesis", ylim=c(12,122))

par(mfrow=c(1,1))

Para la pendiente sacral observamos que las medias de los grupos difieren muy ligeramente, siendo superior en el grupo “Spondylolisthesis”. Existen outliers que pueden alterar la media de los datos, sobretodo en el grupo “Spondylolisthesis”, donde se observan varios puntos que se alejan significativamente del comportamiento del resto de las muestras. Es por eso que probablemente no existan diferencias significativas entre la media de los grupos para esta variable, ya que sin estos outliers la media del grupo “Spondylolisthesis” podría disminuir y se asemejaría aún más a la de los otros grupos.

Variable 5: Radio pélvico (pelvic_radius)

par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Radio pélvico", xlab="Normal", ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_2=="Abnormal"],col="Tomato",main="Radio pélvico", xlab="Abnormal",ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Hernia"],col="coral1",main="Radio pélvico", xlab="Hernia", ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Radio pélvico", xlab="Spondylolisthesis", ylim=c(69,164))

par(mfrow=c(1,1))

Las medias de todos los grupos son semejantes, siendo ligeramente inferiores aquellas del grupo abnormal (incluyendo a los grupos hernia y espondilolistesis por separado). La gran variabilidad en el grupo de espondilolistesis (fíjese en su gran desviación estándar) hace dudar de la significación de esta comparación, aunque sea posible que la variable “radio pélvico” sea distintiva en el grupo “Normal” respecto al resto. Aún así, debido a la variabilidad no esperamos que se cumpla esta hipótesis a priori.

Variable 6: Grado de espondilolistesis (degree_spondylolisthesis)

par(mfrow=c(1,4))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="G.espondilolistesis", xlab="Normal", ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_2=="Abnormal"],col="Tomato",main="G.espondilolistesis", xlab="Abnormal",ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Hernia"],col="coral1",main="G.espondilolistesis", xlab="Hernia", ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="G.espondilolistesis", xlab="Spondylolisthesis", ylim=c(-12,419))

par(mfrow=c(1,1))

Efectivamente el grado de espondilolistesis en pacientes con espondilolistesis es superior respecto al resto, el cual varía ligeramente de la media 0. Esta comparación hace clara la presencia de un elemento de la muestra probablemente mal tomado (existe un paciente con 400 mm (40 cm) de desplazamiento vertebral, hecho bastante cuestionable). No obstante, la comparativa sigue vigente aun ignorar los outliers que distorsionan la media.

Por otro lado, fíjese que esta variable es la única que se mide en milímetros, mostrándonos que la media de variación de la distancia entre vértebras lumbares en pacientes con espondilolistesis es de 51.89 mm. Por lo tanto, la mayoría de pacientes con espondilolistesis del estudio tienen un grado de espondilolistesis positivo, es decir, la vértebra desplazada tiene un desplazamiento anterior respecto su vértebra inferior en la mayoría de casos.

Ejercicio 5:

Para realizar la regresión lineal eligiremos las dos variables que parezcan más relacionadas de entre todas las disponibles. Para ver el comportamiento que cada variable tiene con las otras variables del estudio, utilizaremos un gráfico de dispersión de pares de variables (pairs). Nótese que eliminaremos del análisis las dos variables de tipo caracter (class_2 y class_3) debido a que no son numéricas. Para analizar la relación entre variables de forma cuantitativa, calcularemos las correlaciones entre las variables numéricas mediante el comando “cor()”.

pairs(ortho_df[, -c(7:8)])

cor(ortho_df[, -c(7:8)])
##                          pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## pelvic_incidence                1.0000000  0.62919877            0.71728236
## pelvic_tilt                     0.6291988  1.00000000            0.43276386
## lumbar_lordosis_angle           0.7172824  0.43276386            1.00000000
## sacral_slope                    0.8149600  0.06234529            0.59838689
## pelvic_radius                  -0.2474672  0.03266781           -0.08034361
## degree_spondylolisthesis        0.6387427  0.39786228            0.53366701
##                          sacral_slope pelvic_radius degree_spondylolisthesis
## pelvic_incidence           0.81495999   -0.24746721               0.63874275
## pelvic_tilt                0.06234529    0.03266781               0.39786228
## lumbar_lordosis_angle      0.59838689   -0.08034361               0.53366701
## sacral_slope               1.00000000   -0.34212835               0.52355746
## pelvic_radius             -0.34212835    1.00000000              -0.02606501
## degree_spondylolisthesis   0.52355746   -0.02606501               1.00000000
# Para que sea más facil interpretar las relaciones entre las variables vamos a crear una visualización de las interacciones en contexto de classificación Normal-Anormal y Normal-Hernia-Espondilolistesis.

model.matrix(~0. + pelvic_incidence+pelvic_tilt+pelvic_radius+sacral_slope+lumbar_lordosis_angle+ degree_spondylolisthesis + class_3, ortho_df) %>%
  cor() %>% 
  ggcorrplot(lab=T,
             lab_size=2.5,
             type = "lower")

model.matrix(~0. + pelvic_incidence+pelvic_tilt+pelvic_radius+sacral_slope+lumbar_lordosis_angle+ degree_spondylolisthesis + class_3, ortho_df) %>%
  cor() %>% 
  ggcorrplot(lab=T,
             lab_size=2.5,
             type = "lower")

Este comando (cor()) nos muestra los coeficientes de correlación entre cada variable, por defecto mediante el método de Pearson. Cuando más próximos a 1 o -1, más relación tienen las dos variables, y por lo tanto, el incremento (en caso de coeficiente positivo) o el decremento (en caso de coeficiente negativo) de una variable influencia al incremento/decremento de la otra variable, respectivamente.

Elección de variables De todos los coeficientes de correlación de Pearson entre dos variables distintas, los que tienen la cifra más elevada en valor absoluto son la incidencia pélvica (pelvic_incidence) con la pendiente sacral (sacral_slope), con un coeficiente de 0.8149600; seguido de la incidencia pélvica (pelvic_incidence) con el ángulo de lordosis lumbar (lumbar_lordosis_angle), con con un coeficiente de correlación de 0.71728236.

Hace falta remarcar que existen algunas asociaciones de variables más con coeficientes de correlación elevados que quizás sean importantes en el siguiente análisis de regresión múltiple (véase el Ejercicio 6). No obstante, en este ejercicio, queremos comprobar si es cierto que la incidencia pélvica es una variable que se correlaciona con la pendiente sacral. Al tener un coeficiente de correlación tan próximo a 1, nuestra hipósis es que sí que se cumple esta relación.

De las comparaciones visuales anteriores, si nos fijamos en el gráfico de comparación entre estas dos variables, observamos cómo la nube de puntos de los datos de nuestra muestra forma una línea en diagonal ascendente. Este hecho verifica que el coeficiente de correlación sea tan elevado y que éste tenga un valor positivo: cuando una de las variables aumenta, parece ser que la otra también.

Modelo de regresión lineal Para determinar si esta relación es significativa realizaremos un análisis de regresión lineal entre estas dos variables. Para ello utilizaremos un modelo de regresión lineal (lm) que denominaremos como “ortho_lm”. Nuestra hipótesis será que la pendiente sacrial es la variable independiente (explicativa) y que la incidencia pélvica es la variable dependiente (explicada), ya que el hecho de tener una pendiente sacral más elevada implicaría que el sacro de la persona en cuestión se inclinara hacia atrás y por lo tanto ésta tendría posiblemente un ángulo de incidencia pélvica superior.

ortho_lm <-lm(ortho_df$pelvic_incidence~ortho_df$sacral_slope)
summary(ortho_lm)
## 
## Call:
## lm(formula = ortho_df$pelvic_incidence ~ ortho_df$sacral_slope)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.359  -6.566  -1.403   4.383  32.989 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            15.5461     1.9079   8.148 9.35e-15 ***
## ortho_df$sacral_slope   1.0465     0.0424  24.680  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.01 on 308 degrees of freedom
## Multiple R-squared:  0.6642, Adjusted R-squared:  0.6631 
## F-statistic: 609.1 on 1 and 308 DF,  p-value: < 2.2e-16

Mediante el “summary” del modelo observamos distintos valores estadísticos. El valor principal que nos interesa es el coeficiente de determinación (el R-squared) del modelo, que nos dictará qué tan cerca están los datos de una regresión lineal. En este caso, el coeficiente de determinación no es muy próximo a 1 o -1 (es 0.6642) por lo que la bondad del modelo no parece tan fuerte como creíamos. Según este coeficiente, este modelo explicaría tan solo un 66.42% de la variabilidad de los datos.

Normalidad y varianza de los errores Para comprobar si el coeficiente obtenido es válido, analizaremos la normalidad de los errores del modelo mediante un histograma básico, un boxplot y el gráfico de cuartiles de los residuos.

residuos <-rstandard(ortho_lm)
par(mfrow=c(1,3))
hist(residuos, main = "Histograma de residuos")
boxplot(residuos, main = "Diagrama de cajas de residuos")
qqnorm(residuos, main = "G. cuartiles de residuos")
qqline(residuos)

par(mfrow=c(1,1))

En el primer gráfico observamos cómo, a nivel visual, los residuos parecen tener un comportamiento normal (la forma de una pendiente simétrica con un pico próximo al cero). El segundo gráfico nos añade la información relativa a los outliers. Parece ser que existen diversos valores de nuestra muestra que hacen alejarse los errores de la normalidad de los datos. No obstante, según el tercer gráfico, no parece ser significativo, ya que los residuos siguen una recta diagonal con una desviación en el tercer cuartil únicamente. Así, podemos concluir que los residuos siguen una distribución normal.

Una vez verificada la normalidad de los residuos, analizaremos si la varianza de los errores es constante. El hecho de que la varianza de los errores sea constante implicaría que la variable dependiente variaría debido a nuestra variable independiente, no debido a otros factores externos a nuestro análisis. Para ello, compararemos el gráfico de nuestro modelo con los residuos estandarizados.

plot(fitted.values(ortho_lm),rstandard(ortho_lm),xlab="Valores ajustados",ylab="Residuos estandarizados", col="darkorchid3")
abline(h=0)

Verificamos que existen algunas anotaciones de nuestra muestra (outliers) que se alejan del comportamiento de la mayoría (es decir, que se alejan de la línea centrada al cero). Aun así, éstas parecen no ser significativas debido a que son pocas.

Es por eso que aceptamos el modelo de regresión lineal y suponemos que el hecho de que su coeficiente de determinación no llegue a ser tan elevado como su coeficiente de correlación podría ser debido a la influencia de estos outliers.

Modelo final y conclusiones Finalmente, creemos que el modelo estudiado parece válido y lo representaremos en el siguiente gráfico. Véase cómo la mayoría de datos de nuestra muestra para las dos variables sigue un mismo comportamiento lineal ascendente.

plot(ortho_df$sacral_slope,ortho_df$pelvic_incidence, main="Modelo de regresión lineal", xlab="Pendiente sacral", ylab= "Incidencia pélvica")
abline(ortho_lm)

Como conclusión podemos decir que, aunque el coeficiente de determinación no sea muy elevado, las dos variables analizadas (incidencia pélvica y pendiente sacral) tienen cierta relación en cuanto a su comportamiento. Además, sabemos que esta relación es positiva (un aumento de la pendiente sacral puede implicar un aumento de la incidencia pélvica). Estos datos son consistentes con el comportamiento biológico esperado, ya que el hecho de tener una pendiente sacral más elevada significa que el ángulo entre el sacro y el eje horizontal aumenta, haciendo inclinar así la pelvis. Al hacerlo, el ángulo de incidencia pélvica aumentaría debido a que éste tiene su base de medida entre el eje bicoxofemoral (que no habría cambiado) y línea perpendicular a la placa sacra (que se habría movido).

Este argumento se sustenta de bibliografía, como el estudio de E.V. Geiger et al. (2007), que econtró una correlación entre estas dos variables en su muestra de pacientes.

https://www.researchgate.net/publication/7022399_Adjustment_of_pelvispinal_parameters_preserves_the_constant_gravity_line_position

Ejercicio 6:

Con el análisis del ejercicio 5 hemos podido responder parcialmente la pregunta 2 que nos proponíamos: ¿existe alguna relación entre nuestras variables de estudio? Mediante la regresión lineal hemos comprobado que existe una correlación positiva significativa entre la pendiente sacral y la incidencia pélvica. No obstante, también hemos visto unos índices de correlación elevados entre otras variables de estudio. Es por eso que ahora realizaremos una regresión múltiple para analizar qué variables tienen relación directa con otra/s.

Elección de variables Para definir un modelo inicial de regresión múltiple, debemos indagar de nuevo en los coeficientes de correlación de nuestras variables, para así elegir aquellos que sean próximos a 1 o -1. Nótese que trabajaremos con el mismo método de Pearson para calcularlos.

cor(ortho_df[, -c(7:8)])
##                          pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## pelvic_incidence                1.0000000  0.62919877            0.71728236
## pelvic_tilt                     0.6291988  1.00000000            0.43276386
## lumbar_lordosis_angle           0.7172824  0.43276386            1.00000000
## sacral_slope                    0.8149600  0.06234529            0.59838689
## pelvic_radius                  -0.2474672  0.03266781           -0.08034361
## degree_spondylolisthesis        0.6387427  0.39786228            0.53366701
##                          sacral_slope pelvic_radius degree_spondylolisthesis
## pelvic_incidence           0.81495999   -0.24746721               0.63874275
## pelvic_tilt                0.06234529    0.03266781               0.39786228
## lumbar_lordosis_angle      0.59838689   -0.08034361               0.53366701
## sacral_slope               1.00000000   -0.34212835               0.52355746
## pelvic_radius             -0.34212835    1.00000000              -0.02606501
## degree_spondylolisthesis   0.52355746   -0.02606501               1.00000000

En este caso observamos que, a parte de la correlación significativa entre la incidencia pélvica y la pendiente sacral observada anteriormente, existen asociaciones con coeficientes de correlación elevados.

Es el caso de la incidencia pélvica (pelvic_incidence) con el ángulo de lordosis lumbar (lumbar_lordosis_angle) que tienen un coeficiente de correlación de 0.71728236. Seguidamente, tenemos la incidencia pélvica con el grado de espondilolistesis (degree_spondylolisthesis) (coeficiente de 0.6387427) y la incidencia pélvica con la inclinación pélvica (pelvic_tilt) (coeficiente de 0.6291988). Finalmente existe un coeficiente de correlación de 0.59838689 entre el ángulo de lordosis lumbar y la pendiente sacral (sacral_slope), un coeficiente de 0.53366701 entre el ángulo de lordosis lumbar y el grado de espondilolistesis, y un coeficiente de 0.52355746 entre la pendiente sacral y el grado de espondilolistesis.

A partir de estos casos, tenemos coeficientes de correlación en valor absoluto inferiores a 0.5 por lo que no los incluiremos a priori como destacables.

Por lo tanto, observamos que la variable de incidencia pélvica es la que parece tener más relación con las otras, teniendo los coeficientes de correlación de Pearson más elevados para cuatro de las otras variables. Es por eso que parece intuitivo que esta variable será nuestra variable objetivo para comparar con las otras. En este caso, incluiremos todas las otras variables en el modelo, aunque sabemos que probablemente el radio pélvico sea descartado más adelante, ya que es la única variable de nuestro estudio sin un coeficiente de correlación elevado.

Modelo de regresión múltiple Realizamos el modelo de regresión de la incidencia pélvica con las cinco variables numéricas restantes.

ortho_modelomu<-lm(pelvic_incidence~pelvic_tilt+lumbar_lordosis_angle+sacral_slope+pelvic_radius+degree_spondylolisthesis, data=ortho_df)

summary(ortho_modelomu)
## 
## Call:
## lm(formula = pelvic_incidence ~ pelvic_tilt + lumbar_lordosis_angle + 
##     sacral_slope + pelvic_radius + degree_spondylolisthesis, 
##     data = ortho_df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.065e-08 -4.684e-10 -9.070e-11  4.062e-10  1.100e-08 
## 
## Coefficients:
##                            Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)              -4.423e-10  3.168e-09 -1.400e-01    0.889    
## pelvic_tilt               1.000e+00  3.263e-11  3.065e+10   <2e-16 ***
## lumbar_lordosis_angle     1.887e-11  2.104e-11  8.970e-01    0.371    
## sacral_slope              1.000e+00  3.038e-11  3.292e+10   <2e-16 ***
## pelvic_radius             6.219e-12  2.192e-11  2.840e-01    0.777    
## degree_spondylolisthesis -5.313e-12  9.453e-12 -5.620e-01    0.575    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.688e-09 on 304 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 8.354e+20 on 5 and 304 DF,  p-value: < 2.2e-16

El “summary” del modelo nos informa de la significación de la correlación de nuestra variable objetivo (incidencia pélvica) con cada otra variable. Obsérvese que los coeficientes significativos son únicamente con la inclinación pélvica (pelvic_tilt), con p-valor inferior a 2e-16, y con la pendiente sacral (sacral_slope), con p-valor inferior a 2e-16. Los asteriscos nos informan visualmente de que éstas son las únicas dos variables con correlaciones significativas para la incidencia pélvica.

Además, nos otorga el coeficiente de determinación (R-squared) del modelo, siendo éste el máximo (1). Esta información parece sospechosa, pues en biología es difícil obtener un coeficiente de determinación total, pero si nos fijamos en el p-valor del modelo, observamos que éste es significativo (inferior a 2.2e-16). Es probable que tengamos un caso de multicolinealidad, es decir, varias variables que tienen correlación muy elevada, y eso haría incrementar el coeficiente de determinación. Aun así, como es un análisis significativo, ésto nos verifica que la bondad del modelo de regresión es positiva. En otras palabras, toda la variabilidad de la variable en cuestión (incidencia pélvica) viene dada por otras variables del estudio.

Ahora que sabemos que nuestras variables predictoras para la incidencia pélvica son la pendiente sacral y la inclinación pélvica, perfeccionaremos el modelo de regresión.

ortho_modelomu2<-lm(pelvic_incidence~pelvic_tilt+sacral_slope, data=ortho_df)
summary(ortho_modelomu2)
## 
## Call:
## lm(formula = pelvic_incidence ~ pelvic_tilt + sacral_slope, data = ortho_df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.038e-08 -3.259e-10 -7.410e-11  4.869e-10  1.077e-08 
## 
## Coefficients:
##               Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)  7.766e-10  9.827e-10 7.900e-01     0.43    
## pelvic_tilt  1.000e+00  2.662e-11 3.757e+10   <2e-16 ***
## sacral_slope 1.000e+00  1.985e-11 5.039e+10   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.674e-09 on 307 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.101e+21 on 2 and 307 DF,  p-value: < 2.2e-16

En este sentido, si aceptamos el modelo, la variable “incidencia pélvica” tendría toda su variabilidad explicada por la pendiente sacral (como hemos corroborado en el ejercicio anterior) y la inclinación pélvica.

A nivel anatómico, estos datos tienen mucho sentido. Como se ha visto en el ejercicio anterior, existe una relación documentada entre la incidencia pélvica y la pendiente sacral. En cuanto a la incidencia pélvica y la inclinación pélvica, pasa exactamente lo mismo.

De hecho, tal como indica el artículo de J.C.Le Huec et al. (2011), en medicina la incidencia pélvica (PI) se calcula mediante la suma de la pendiente sacral (SS) más la inclinación pélvica (PT).

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3175921/

Por lo tanto, parece lógico que estas tres variables tengan una correlación positiva significativa y que, por lo tanto, a mayor inclinación pélvica, mayor sea la incidencia pélvica.

Ejerecicio 7:

Para decidir si se puede utuilizar el test ANOVA, tenemos que confirmar inialmente si las variables sigue la distribución normal y si la variable no es homoscedastica (tiene la misma varianza).

# Empezamos mmirando a las variables y calculando valores medias y error estandart. 
ortho_df2 %>% group_by(type_of_measure) %>% summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))
## # A tibble: 5 × 3
##   type_of_measure       `Mean of angle` `Standart deviation`
##   <chr>                           <dbl>                <dbl>
## 1 lumbar_lordosis_angle            51.9                 18.6
## 2 pelvic_incidence                 60.5                 17.2
## 3 pelvic_radius                   118.                  13.3
## 4 pelvic_tilt                      17.5                 10.0
## 5 sacral_slope                     43.0                 13.4
ortho_df2 %>%  summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))
## # A tibble: 1 × 2
##   `Mean of displacement` `SD of displacement`
##                    <dbl>                <dbl>
## 1                   26.3                 37.5
# Creamos la visualización de la distribución de las variables con denisity plot y con boxplot, marcando valores médias. En la ditribución normal valor media y valor mediana y moda coinciden. Por la visualización, parece que la variable pelvic_radius y la sacral_slope pueden seguir la distribución normal. También miramos si el hecho de pertenecer a una otra classificación puede afectar drasticamente al tipo de la sistribución. Se evidencia el cambio importante de comportamiento de la variable degree_spondilolistesis en comparación normal/anormal y espondilolistesis vs el resto de classificaciones.

pt <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)
ptc2 <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)+ facet_grid(~class_2)
ptc3 <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)+ facet_grid(~class_3)
pp <- list(pt,ptc2,ptc3)
ggarrange(nrow = 3, ncol= 1, plotlist = pp)

pds <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis)) + geom_density(alpha = 0.2, fill = "darkgreen")
pds2 <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis, fill = class_2)) + geom_density(alpha = 0.2)+ facet_grid(~class_2)
pds3 <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis, fill = class_3)) + geom_density(alpha = 0.2)+ facet_grid(~class_3)
ppds <- list(pds,pds2,pds3)
ggarrange(nrow = 3, ncol= 1, plotlist = ppds)

# Evaluamos las valores medias y del error estandar segun la classificación:
# Segun class_2
ortho_df2 %>% group_by(type_of_measure, class_2) %>%
  summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))
## `summarise()` has grouped output by 'type_of_measure'. You can override using the `.groups` argument.
## # A tibble: 10 × 4
## # Groups:   type_of_measure [5]
##    type_of_measure       class_2  `Mean of angle` `Standart deviation`
##    <chr>                 <fct>              <dbl>                <dbl>
##  1 lumbar_lordosis_angle Normal              43.5                12.4 
##  2 lumbar_lordosis_angle Abnormal            55.9                19.7 
##  3 pelvic_incidence      Normal              51.7                12.4 
##  4 pelvic_incidence      Abnormal            64.7                17.7 
##  5 pelvic_radius         Normal             124.                  9.01
##  6 pelvic_radius         Abnormal           115.                 14.1 
##  7 pelvic_tilt           Normal              12.8                 6.78
##  8 pelvic_tilt           Abnormal            19.8                10.5 
##  9 sacral_slope          Normal              38.9                 9.62
## 10 sacral_slope          Abnormal            44.9                14.5
ortho_df2 %>% group_by(class_2) %>%
  summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))
## # A tibble: 2 × 3
##   class_2  `Mean of displacement` `SD of displacement`
##   <fct>                     <dbl>                <dbl>
## 1 Normal                     2.19                 6.28
## 2 Abnormal                  37.8                 40.6
# Según class_3
ortho_df2 %>% group_by(type_of_measure, class_3) %>% 
  summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))
## `summarise()` has grouped output by 'type_of_measure'. You can override using the `.groups` argument.
## # A tibble: 15 × 4
## # Groups:   type_of_measure [5]
##    type_of_measure       class_3           `Mean of angle` `Standart deviation`
##    <chr>                 <fct>                       <dbl>                <dbl>
##  1 lumbar_lordosis_angle Normal                       43.5                12.4 
##  2 lumbar_lordosis_angle Hernia                       35.5                 9.77
##  3 lumbar_lordosis_angle Spondylolisthesis            64.1                16.4 
##  4 pelvic_incidence      Normal                       51.7                12.4 
##  5 pelvic_incidence      Hernia                       47.6                10.7 
##  6 pelvic_incidence      Spondylolisthesis            71.5                15.1 
##  7 pelvic_radius         Normal                      124.                  9.01
##  8 pelvic_radius         Hernia                      116.                  9.36
##  9 pelvic_radius         Spondylolisthesis           115.                 15.6 
## 10 pelvic_tilt           Normal                       12.8                 6.78
## 11 pelvic_tilt           Hernia                       17.4                 7.02
## 12 pelvic_tilt           Spondylolisthesis            20.7                11.5 
## 13 sacral_slope          Normal                       38.9                 9.62
## 14 sacral_slope          Hernia                       30.2                 7.56
## 15 sacral_slope          Spondylolisthesis            50.8                12.3
ortho_df2 %>% group_by(class_3) %>%  
  summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))
## # A tibble: 3 × 3
##   class_3           `Mean of displacement` `SD of displacement`
##   <fct>                              <dbl>                <dbl>
## 1 Normal                              2.19                 6.28
## 2 Hernia                              2.48                 5.49
## 3 Spondylolisthesis                  51.9                 40.0
# Creamos la visualización de valores médias y error estandar de las variables.
angles_ms <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
  stat_summary(fun = mean,
               shape = 18,
               size = 1)+ 
  stat_summary(fun = mean,
               fun.min = function(x) {mean(x) - sd(x)},
               fun.max =function(x) {mean(x) + sd(x)},
               geom = "errorbar")+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
  
angles_ms2 <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
  stat_summary(fun = mean,
               shape = 18,
               size = 1)+ 
  stat_summary(fun = mean,
               fun.min = function(x) {mean(x) - sd(x)},
               fun.max =function(x) {mean(x) + sd(x)},
               geom = "errorbar")+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank()) +
  facet_wrap(~class_2)

angles_ms3 <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
  stat_summary(fun = mean,
               shape = 18,
               size = 1)+ 
  stat_summary(fun = mean,
               fun.min = function(x) {mean(x) - sd(x)},
               fun.max =function(x) {mean(x) + sd(x)},
               geom = "errorbar")+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank()) +
  facet_wrap(~class_3)

degree_ms <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = 0)) +
  stat_summary(fun = mean,
               shape = 18,
               size = 1,
               color = "green")+ 
  stat_summary(fun = mean,
               fun.min = function(x) {mean(x) - sd(x)},
               fun.max =function(x) {mean(x) + sd(x)},
               color = "green",
               geom = "errorbar")+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())

degree_ms2 <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = class_2, color = class_2)) +
  stat_summary(fun = mean,
               shape = 18,
               size = 1)+ 
  stat_summary(fun = mean,
               fun.min = function(x) {mean(x) - sd(x)},
               fun.max =function(x) {mean(x) + sd(x)},
               geom = "errorbar")+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())  +
  facet_wrap(~class_2)

degree_ms3 <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = class_3, color = class_3)) +
  stat_summary(fun = mean,
               shape = 18,
               size = 1)+ 
  stat_summary(fun = mean,
               fun.min = function(x) {mean(x) - sd(x)},
               fun.max =function(x) {mean(x) + sd(x)},
               geom = "errorbar")+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank()) +
  facet_wrap(~class_3)

ggarrange(plotlist= list(angles_ms, angles_ms2, angles_ms3) , nrow= 3, ncol = 1, legend = "left")
## Warning: Removed 5 rows containing missing values (geom_segment).

## Warning: Removed 5 rows containing missing values (geom_segment).

## Warning: Removed 5 rows containing missing values (geom_segment).

## Warning: Removed 5 rows containing missing values (geom_segment).

## Warning: Removed 5 rows containing missing values (geom_segment).

## Warning: Removed 5 rows containing missing values (geom_segment).

ggarrange(plotlist= list(degree_ms, degree_ms2, degree_ms3) , nrow= 3, ncol = 1, legend = "left")
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).

## Warning: Removed 1 rows containing missing values (geom_segment).

## Warning: Removed 1 rows containing missing values (geom_segment).

## Warning: Removed 1 rows containing missing values (geom_segment).

## Warning: Removed 1 rows containing missing values (geom_segment).

#Resumimos la data en el grafico tipo boxplot para evaluar los quartiles y su relación con valores médias:

ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, fill = type_of_measure)) +
  geom_boxplot() + stat_summary(fun = mean,
                                shape = 18,
                                size = 1) + 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
## Warning: Removed 5 rows containing missing values (geom_segment).

ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure), color = "black") +
  geom_boxplot(alpha = 0.3, aes(fill = type_of_measure)) +
  stat_summary(fun.y = mean,
               geom = "point",
               shape = 18, size = 4,
               aes(color = type_of_measure))+
  facet_wrap(~class_2)+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
## Warning: `fun.y` is deprecated. Use `fun` instead.

ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure), color = "black") +
  geom_boxplot(alpha = 0.3, aes(fill = type_of_measure)) +
  stat_summary(fun.y = mean, 
               geom = "point", shape = 18,
               size = 4, 
               aes(color = type_of_measure)) +
  facet_wrap(~class_3)+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
## Warning: `fun.y` is deprecated. Use `fun` instead.

ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis,
                         x = 0)) +
  geom_boxplot(fill = "red",
               alpha = 0.4) +
  stat_summary(fun = mean, shape = 18,
               size =1)+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
## Warning: Removed 1 rows containing missing values (geom_segment).

ortho_df2 %>% ggplot(aes(x =0,
                         y= degree_spondylolisthesis,
                         fill= class_2),
                     color = black) +
  geom_boxplot(alpha = 0.6,
               aes(fill = class_2)) +
  stat_summary(fun.y = mean,
               geom = "point",
               shape = 18, size = 5, 
               aes(color = class_2))+
  facet_wrap(~class_2)+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
## Warning: `fun.y` is deprecated. Use `fun` instead.

ortho_df2 %>% ggplot(aes(x = 0, y= degree_spondylolisthesis), color = black) +
  geom_boxplot(alpha = 0.4, aes(fill = class_3)) +
  stat_summary(fun.y = mean, geom = "point", shape = 18, size = 5, aes(color = class_3))+ facet_wrap(~class_3)+ 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank())
## Warning: `fun.y` is deprecated. Use `fun` instead.

# Normality tests: Ahora vamos a comprobar si la distribución de las variables realmente sigue la distribución normal (hipotesis 0). Para evaluarlo vamos a utilizar el test de Kolmogorov-Smirnov, y la visualización con QQ-plot

# Kolmogorov-Smirnov
# Para realizar este test sobre el conjunto de data vamos a utilizar función map - una varaante de lapply() dentro del paquete tidyverse (libreria purr)

ortho_df %>% 
  select(where(is.numeric))%>% 
  colnames() %>%
  set_names() %>% 
  map(~ lillie.test(ortho_df[,.])) %>%
  map_dfr(., tidy, .id = "variable")
## # A tibble: 6 × 4
##   variable                 statistic  p.value method                            
##   <chr>                        <dbl>    <dbl> <chr>                             
## 1 pelvic_incidence            0.0708 7.33e- 4 Lilliefors (Kolmogorov-Smirnov) n…
## 2 pelvic_tilt                 0.0829 2.37e- 5 Lilliefors (Kolmogorov-Smirnov) n…
## 3 lumbar_lordosis_angle       0.0687 1.25e- 3 Lilliefors (Kolmogorov-Smirnov) n…
## 4 sacral_slope                0.0578 1.43e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 5 pelvic_radius               0.0593 1.06e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 6 degree_spondylolisthesis    0.173  5.13e-25 Lilliefors (Kolmogorov-Smirnov) n…
# Vamos a mirar si el hecho de quitar el outlier (valor 418 del grado de espondilolistesis nos afecta al tipo de distribución)
ortho_df3 %>% 
  select(where(is.numeric))%>% 
  colnames() %>%
  set_names() %>% 
  map(~ lillie.test(ortho_df[,.])) %>%
  map_dfr(., tidy, .id = "variable")  # Parece que no afecta, así que procedemos a seguir con el análisis con el data.frame completo.
## # A tibble: 6 × 4
##   variable                 statistic  p.value method                            
##   <chr>                        <dbl>    <dbl> <chr>                             
## 1 pelvic_incidence            0.0708 7.33e- 4 Lilliefors (Kolmogorov-Smirnov) n…
## 2 pelvic_tilt                 0.0829 2.37e- 5 Lilliefors (Kolmogorov-Smirnov) n…
## 3 lumbar_lordosis_angle       0.0687 1.25e- 3 Lilliefors (Kolmogorov-Smirnov) n…
## 4 sacral_slope                0.0578 1.43e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 5 pelvic_radius               0.0593 1.06e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 6 degree_spondylolisthesis    0.173  5.13e-25 Lilliefors (Kolmogorov-Smirnov) n…
### QQplot

ggqqplot(ortho_df2, "angles", facet.by = "class_3", color = "type_of_measure")

ggqqplot(ortho_df2, "angles", facet.by = "class_2", color = "type_of_measure")

ggqqplot(ortho_df2, "degree_spondylolisthesis", color = "class_2", facet.by  = "class_2")

ggqqplot(ortho_df2, "degree_spondylolisthesis", color = "class_3", facet.by  = "class_3")

### Homogenidad - vamos a mirar si la data es homogénea
ortho_df %>% 
  select(where(is.numeric))%>% 
  colnames() %>%
  set_names() %>%  
  map(~ fligner.test(ortho_df[,.] ~ class_2, data = ortho_df)) %>%
  map_dfr(., tidy, .id = "variable")
## # A tibble: 6 × 5
##   variable                 statistic  p.value parameter method                  
##   <chr>                        <dbl>    <dbl>     <dbl> <chr>                   
## 1 pelvic_incidence             18.9  1.40e- 5         1 Fligner-Killeen test of…
## 2 pelvic_tilt                  14.0  1.80e- 4         1 Fligner-Killeen test of…
## 3 lumbar_lordosis_angle        25.3  4.98e- 7         1 Fligner-Killeen test of…
## 4 sacral_slope                 16.4  5.11e- 5         1 Fligner-Killeen test of…
## 5 pelvic_radius                 9.29 2.30e- 3         1 Fligner-Killeen test of…
## 6 degree_spondylolisthesis     89.8  2.58e-21         1 Fligner-Killeen test of…
ortho_df %>% 
  select(where(is.numeric))%>% 
  colnames() %>%
  set_names() %>%  
  map(~ fligner.test(ortho_df[,.] ~ class_3, data = ortho_df)) %>%
  map_dfr(., tidy, .id = "variable")
## # A tibble: 6 × 5
##   variable                 statistic  p.value parameter method                  
##   <chr>                        <dbl>    <dbl>     <dbl> <chr>                   
## 1 pelvic_incidence              7.21 2.72e- 2         2 Fligner-Killeen test of…
## 2 pelvic_tilt                  27.3  1.19e- 6         2 Fligner-Killeen test of…
## 3 lumbar_lordosis_angle        19.3  6.56e- 5         2 Fligner-Killeen test of…
## 4 sacral_slope                  7.28 2.63e- 2         2 Fligner-Killeen test of…
## 5 pelvic_radius                22.9  1.07e- 5         2 Fligner-Killeen test of…
## 6 degree_spondylolisthesis    125.   6.04e-28         2 Fligner-Killeen test of…

Anova test:

Tras realizar varios tests de normalidad se confirma que la distribución de data NO es normal, asi que procedemos a realizar tests ANOVA, con la idéa en mente de que la distribución puede ser insuficientemente normal para conseguir resultados fiables. El test de homogenidad de data sugiere NO-homogenidad de data.

# Realizoms el test de anova para diferentes variables.
anova_PI <- aov(pelvic_incidence~class_3, ortho_df)
anova_S <- aov(degree_spondylolisthesis~class_3, ortho_df)
anova_PR <- aov(pelvic_radius~class_3, ortho_df)
anova_LLA <- aov(lumbar_lordosis_angle~class_3, ortho_df)
# La relación de las medias entre gurpos para la incidencia pélvica
par(mfrow=c(2,2))
plot(anova_PI)

par(mfrow=c(1,1))
plot(TukeyHSD(anova_PI))

# La relación de las medias entre gurpos para el grado de espondilolistesis
par(mfrow=c(2,2))
plot(anova_S)

par(mfrow=c(1,1))
plot(TukeyHSD(anova_S))

# La relación de las medias entre gurpos para angulo del radius de pelvis
par(mfrow=c(2,2))
plot(anova_PR)

par(mfrow=c(1,1))
plot(TukeyHSD(anova_PR))

# La relación de las medias entre gurpos para el angulo del lordosis lumbar
par(mfrow=c(2,2))
plot(anova_LLA)

par(mfrow=c(1,1))
plot(TukeyHSD(anova_LLA))

Los tests sugieren que la diferencia de los valores medias de algunas variables se difiere de forma significativa entre algunos grupos. (ver ejercicio 8)

 #evaluamos un scatter plot para evaluar presencia de algun tipo de aglomeración.
ortho_df %>% ggplot(aes(lumbar_lordosis_angle, pelvic_incidence, color = class_3)) +
  geom_point()

ortho_df %>% ggplot(aes(lumbar_lordosis_angle, pelvic_incidence, color = class_2)) +
  geom_point()

# Miramos si se puede identificar algún numero de cluster
fviz_nbclust(ortho_df[,c(1,3)], kmeans, method = "wss")

fviz_nbclust(ortho_df[,c(1,3)], kmeans, method = "silhouette")

# Probamos classificar las valores segun la distancia entre los puntos del scatter plot para entre 2 y 3 clusteres

cluster2 <- kmeans(ortho_df[, c(1,3)], 2, nstart = 25)
table(cluster2$cluster, ortho_df$class_2)
##    
##     Normal Abnormal
##   1     16      116
##   2     84       94
summary(cluster2)
##              Length Class  Mode   
## cluster      310    -none- numeric
## centers        4    -none- numeric
## totss          1    -none- numeric
## withinss       2    -none- numeric
## tot.withinss   1    -none- numeric
## betweenss      1    -none- numeric
## size           2    -none- numeric
## iter           1    -none- numeric
## ifault         1    -none- numeric
fviz_cluster(cluster2, data = ortho_df[, c(1,3)])

cluster2 # Se evidencia que la classificación entre 2 clusteres es bastente complicada y la busqueda de algomeraciones por distancia no ha sido muy eficaz, y se acierta solo en 61.4% de casos.
## K-means clustering with 2 clusters of sizes 132, 178
## 
## Cluster means:
##   pelvic_incidence lumbar_lordosis_angle
## 1         76.19685              68.75039
## 2         48.85381              39.45807
## 
## Clustering vector:
##   [1] 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2
## [149] 2 2 2 2 1 2 2 1 2 2 2 1 2 2 2 1 1 1 1 2 2 2 2 2 1 1 2 2 1 1 1 2 2 1 1 1 1
## [186] 2 1 2 1 2 1 1 2 1 2 1 1 2 2 2 2 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 43253.08 33224.06
##  (between_SS / total_SS =  61.4 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
cluster3 <- kmeans(ortho_df[, c(1,3)], 3, nstart = 25)
table(cluster3$cluster, ortho_df$class_3)
##    
##     Normal Hernia Spondylolisthesis
##   1     38     12                67
##   2     58     48                13
##   3      4      0                70
summary(cluster3)
##              Length Class  Mode   
## cluster      310    -none- numeric
## centers        6    -none- numeric
## totss          1    -none- numeric
## withinss       3    -none- numeric
## tot.withinss   1    -none- numeric
## betweenss      1    -none- numeric
## size           3    -none- numeric
## iter           1    -none- numeric
## ifault         1    -none- numeric
fviz_cluster(cluster3, data = ortho_df[, c(1,3)])

cluster3 # y la clasterización en 3 clusteres se acierta en 74.1% de casos.
## K-means clustering with 3 clusters of sizes 117, 119, 74
## 
## Cluster means:
##   pelvic_incidence lumbar_lordosis_angle
## 1         63.58681              53.46109
## 2         43.99372              35.12170
## 3         82.14937              76.54268
## 
## Clustering vector:
##   [1] 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 1 1 2 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 1 2 2 2 1 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2 1 2 2 1 2 1 1 1 2 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 3 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 1 3 3 1 1 1 3 3
## [223] 1 1 1 1 3 1 3 3 1 3 1 1 3 1 1 1 3 3 1 3 1 1 3 1 1 1 3 3 3 3 1 3 3 1 3 3 3
## [260] 3 3 1 1 3 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 1
## [297] 3 3 1 3 3 3 3 3 3 3 3 3 3 3
## 
## Within cluster sum of squares by cluster:
## [1] 17597.78 12236.25 21400.69
##  (between_SS / total_SS =  74.1 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
# otra forma de realizar analisis de cluster en estos casos es clusterización jerárquica
cluster <- hclust(dist(ortho_df), method="complete")
## Warning in dist(ortho_df): NAs introduced by coercion
plot(cluster)

pltree(agnes(ortho_df, method = "complete"))

Ejercicio 8:

Mediante este estudio hemos podido resolver las dos preguntas principales hechas en el principio del trabajo.

Pregunta 1: ¿Existen diferencias significativas entre las variables angulares del grupo “Normal” respecto al grupo “Abnormal”?.

Según la muestra analizada, sí que existe una diferencia significativa entre algunas de las variables según los distintos grupos. Primeramente hemos visto cómo las medias de los pacientes con espondilolistesis para la incidencia pélvica y el ángulo de lordosis lumbar son superiores que en sujetos del grupo “Normal”. Este hecho denota que el desliz de las vértebras de los pacientes con espondilolistesis afecta a estas dos variables métricas más que a los otros sujetos. Por lo tanto, los resultados sugieren que la incidencia pélvica y el ángulo de lordosis lumbar podrían llegar a ser dos variables definitorias para la espondilolistesis. Para los pacientes con “Hernia”, no obstante, no había diferencias significativas respecto a los sujetos del grupo “Normal”. Por lo tanto, el desplazamiento vertebral (espondilolistesis) influye claramente más que la hernia en las variables métricas posturales.

Por otro lado, tenemos que no parecen existir diferencias significativas entre los grupos para el radio pélvico, la pendiente sacral y la inclinación pélvica, por lo que éstas podrían ser variables menos explicativas para las afectaciones estudiadas. De hecho, el grupo de pacientes con Hernia no difiere significativamente del grupo Normal para ninguna variable, dando a entender que las variables métricas estudiadas en este análisis no serían buenos distintivos para esta patología. Quizás con más estudios de muestras mayores y más representativas se podría contrastar los datos de las dos patologías y verificar si las variables utilizadas en estos ejercicios son buenas variables tanto definidoras como métricas.

Pregunta 2: ¿Existe alguna relación o correlación entre las variables del estudio?

Según los resultados de nuestros análisis, sí que existe una relación entre algunas de las variables.

Primeramente, hemos detectado una correlación significativa entre la incidencia pélvica con la pendiente sacral. Esta correlación se ha detectado mediante un análisis de regresión lineal y ha sido confirmado con un análisis de regresión múltiple. Con éste último también se ha detectado una correlación elevada entre la incidencia pélvica y la inclinación pélvica. La relación entre estas tres variables ha sido verificada con los datos en la bibliografía. Además, esta relación tiene sentido a nivel anatómico, ya que las modificaciones en la inclinación pélvica y/o la pendiente sacral hacen que el ángulo de incidencia pélvica cambia cuando se progresa la enfermedad (presencia de hernias o espondilolistesis).

Durante el análisis de los valores médias entre los grupos (ANOVA) se evidencia que, aunque la distribución de las variables no es normal, lo que hace que el test elegido no es muy fiable, en algunos grupos el valor média es diferente de forma significativa.

Variable Relacion
pelvic_incidence NO varia entre el grupo normal y la hernia, pero varia en comparación con el grupo espondilolistesis
degree_spondylolisthesis NO varia entre el grupo normal y la hernia, pero varia en comparación con el grupo espondilolistesis
pelvic_radius NO varia entre el grupo espondilolistesis y la hernia, pero varia en comparación con el grupo espondilolistesis y el grupo normal
lumbar_lordosis_angle Varia de forma importante entre los tres grupos.

Durante el analisi de cluster se realiza 2 tipos de análisis - no jerárquico utilizando valores de distancia entre puntos, y jerárquica. Al realizar el test kmeans y al evaluar el data.frame se evidencia que el numero optimo de clusters se situa entre 2 y 3, por lo que se realiza analisis para las dos situaciones, objetivando mejor calidad de aciertamiento (74%) de clusterización para la situación con 3 clusteres.